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Abstract 


In  this  paper  a  fast  algorithm  for  computing  the  capacitance  of  a  complicated  3-D  geometry 
of  ideal  conductors  in  a  uniform  dielectric  is  described.  The  method  is  an  acceleration  of 
the  standard  integral  equation  approach  for  multiconductor  capacitance  extraction.  These 
integral  equation  methods  are  slow  because  they  lead  to  dense  matrix  problems  which  are 
typically  solved  with  some  form  of  Gaussian  elimination.  This  implies  the  computation 
grows  like  n3,  where  n  is  the  number  of  tiles  needed  to  accurately  discretize  the  conductor 
surface  charges.  In  this  paper  we  present  a  preconditioned  conjugate-gradient  iterative 
algorithm  with  a  multipole  approximation  to  compute  the  iterates.  This  reduces  the 
complexity  so  that  accurate  multiconductor  capacitance  calculations  grow  as  nm  where  m  is 
the  number  of  conductors. 
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Abstract 

In  this  paper  a  fast  algorithm  for  computing  the  capacitance  of  a 
complicated  3-D  geometry  of  ideal  conductors  in  a  uniform  dielectric 
is  described.  The  method  is  an  acceleration  of  the  standard  integral 
equation  approach  for  multiconductor  capacitance  extraction.  These 
integral  equation  methods  are  slow  because  they  lead  to  dense  matrix 
problems  which  are  typically  solved  with  some  form  of  Gaussian  elim¬ 
ination.  This  implies  the  computation  grows  like  n3,  where  n  is  the 
number  of  tiles  needed  to  accuracy  discretize  the  conductor  surface 
charges.  In  this  paper  we  present  a  preconditioned  conjugate-gradient 
iterative  algorithm  with  a  multipole  approximation  to  compute  the 
iterates.  This  reduces  the  complexity  so  that  accurate  multiconduc¬ 
tor  capacitance  calculations  grow  as  nm  where  m  is  the  number  of 
conductors. 


1  Introduction 


In  the  design  of  high  performance  integrated  circuits,  there  are  many  cases 
where  accurate  estimates  of  the  capacitances  of  complicated  three  dimen¬ 
sional  structures  are  important  for  determining  final  circuit  speeds  or  func¬ 
tionality.  Two  examples  are  complicated  three-dimensional  dynamic  mem¬ 
ory  cells  and  the  three-dimensional  chip  carriers  commonly  used  in  main¬ 
frame  computers.  In  these  problems,  capacitance  extraction  is  made  tractable 
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by  assuming  the  conductors  are  ideal,  and  are  embedded  in  a  piecewise- 
constant  dielectric  medium.  Then  to  compute  the  capacitances,  Laplace's 
equation  is  solved  numerically  over  the  charge  free  region  with  the  conduc¬ 
tors  providing  boundary  conditions. 

Although  there  are  a  variety  of  numerical  methods  that  can  be  used  to 
solve  Laplace’s  equation,  the  technique  that  is  typically  used  in  three  di¬ 
mensions  is  the  integral  equation  approach[ruehli73,rao84,ning88].  In  this 
approach,  the  surfaces  or  edges  of  all  the  conductors  are  broken  into  small 
tiles.  It  is  assumed  that  on  each  tile  t,  a  charge,  is  uniformly  or  linearly 
distributed.  The  potential  on  each  tile  is  then  computed  by  summing  the 
contributions  to  the  potential  from  all  the  tiles  using  Laplace’s  equation 
Green’s  functions.  In  this  way  a  matrix  of  potential  coefficients,  P,  relating 
the  set  of  n  tile  potentials  and  the  set  of  n  tile  charges  is  constructed,  and 
must  be  solved  to  compute  capacitances.  Typically,  Guassian  elimination 
or  Cholesky  factorization  is  used  to  solve  the  equation,  in  which  case  the 
number  of  operations  is  order  n3.  Clearly,  this  approach  becomes  compu¬ 
tationally  intractable  if  the  number  of  unknowns  exceeds  several  hundred, 
and  this  limits  the  size  of  the  problem  that  can  be  analyzed  to  one  with  a 
few  conductors. 

In  this  paper  we  present  an  algorithm  for  computing  capacitance  whose 
complexity  grows  as  mn,  where  m  is  the  number  of  conductors.  Our  algo¬ 
rithm,  which  is  really  the  pasting  together  of  three  well-known  algorithms 
[rohklin86],  is  presented  in  three  sections.  To  begin,  in  the  next  section  one 
of  the  standard  integral  equation  approaches  is  briefly  described,  and  it  is 
shown  that  the  algorithm  requires  the  solution  of  an  n  x  n  dense  symmetric 
matrix.  Then,  in  Section  3,  a  preconditioned  conjugate-gradient  algorithm 
is  described,  and  it  is  shown  to  reduce  the  complexity  of  the  calculation  to 
order  mn2.  In  Section  4,  it  is  shown  that  the  conjugate-gradient  algorithm 
only  requires  the  evaluation  of  a  potential  field  from  a  charge  distribution, 
and  this  can  be  computed  in  order  n  time  using  a  multipole  algorithm.  In 
Section  5,  some  preliminary  experimental  results  are  given,  and  we  present 
our  conclusions  and  acknowledgments. 

2  The  Integral  Equation  Approach 

Consider  a  system  of  m  ideal  conductors  embedded  in  a  uniform  lossless 
dielectric  medium.  For  such  a  system,  the  relation  between  the  m  conductor 
potentials,  denoted  by  p  G  and  the  m  total  charges  on  each  conductor, 
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denoted  by  q  €  3?m,  is  given  by  q  =  Cp,  where  C  €  §?mxm  jg  referred  to  as 
the  capacitance  matrix.  The  ith  column  of  C  can  be  calculated  by  solving 
for  the  total  charges  on  each  of  the  conductors  when  the  ith  conductor  is  at 
unit  potential,  and  all  the  other  conductors  are  at  zero  potential.  Then  the 
charge  on  conductor  j,  qj,  is  equal  to  Cij. 

There  are  a  variety  of  approaches  for  numerically  computing  the  con¬ 
ductor  charges  given  a  set  of  conductor  potentials,  and  we  will  focus  on 
integral  equation  methods[ruehli73,rao84,ning88],  as  they  are  efficient  when 
applied  to  problems  with  ideal  conductors  in  a  uniform  dielectric  medium. 
The  method  exploits  the  fact  that  the  charge  is  restricted  to  the  surface  of 
the  conductors,  and  rather  than  discretizing  all  of  free  space,  just  the  sur¬ 
face  charge  on  the  conductors  is  discretized.  The  potential  is  related  to  the 
discretized  surface  charge  through  integrals  of  a  Green’s  functions. 

Let  the  surfaces  of  a  collection  of  m  conductors  in  free  space  be  dis¬ 
cretized  into  a  total  of  n  tiles.  The  potential  at  the  center  of  the  ith  tile 
would  be  the  sum  of  the  contributions  to  the  potential  from  the  charge 
distribution  on  every  tile.  That  is, 


Pi 


(1) 


where  fi  is  the  position  of  the  center  of  tile  i,  r  is  the  position  on  the  surface 
of  tile  j ,  p,  is  the  potential  at  fi,  qj(r)  is  the  position  dependent  charge 
density  on  the  surface  of  the  jth  tile,  and  |  r  |  denotes  the  Euclidian  length  of 
t.  Note  that  the  integral  in  (1)  is  the  free  space  Green’s  function  multiplied 
by  the  charge  density,  integrated  over  the  surface  of  the  jtfl  tile,  and  that  as 
the  distance  between  tile  i  and  tile  j  becomes  large  compared  to  the  surface 
area  of  tile  j,  the  integral  reduces  to  p-tfrj  wbere  qj  is  the  total  charge  on 
tile  j. 

There  are  several  approaches  to  simplifying  (1),  the  simplest  is  the 
“point-matching”  approximation  in  which  it  is  assumed  that  the  charge 
is  distributed  uniformly  on  the  tile  surface[rao84].  In  that  case  (1)  can  be 
simplified  to 


(2) 


where  qj  is  the  total  charge  on  tile  j ,  and  aj  is  the  surface  area  of  tile  j. 
When  applied  to  the  collection  of  n  tiles,  a  dense  linear  system  results, 


Pq  =  p 


(3) 
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where  P  6  §?"xn;  q,  p  €  §?"  and 


Pij  =  Pa  = 

\\-j  pAr-rda+I/  -.-da 
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(4) 


Note  that  q  and  p  are  the  vectors  of  tile  charges  and  potentials  rather  than 
the  conductor  charge  and  potential  vectors,  q  and  p  mentioned  above.  In  (4), 
the  potential  coefficients,  P*  have  been  “symmetrized”  by  averaging  for  sev¬ 
eral  reasons:  the  physical  system  is  symmetric,  the  symmetrized  equations 
have  been  shown  to  produce  more  accurate  results  for  a  given  discretization, 
and  a  symmetric  matrix  problem  is  more  easily  solved.  The  dense  linear 
system  of  (3)  can  be  solved,  typically  by  some  form  of  symmetric  Gaussian 
elimination,  to  compute  tile  charges  from  a  given  set  of  tile  potentials.  To 
compute  the  jth  column  of  the  capacitance  matrix,  (3)  must  be  solved  for 
q ,  given  a  p  vector  whose  entries  p,  are  set  equal  to  one  if  tile  t  is  on  the 
jth  conductor,  and  zero  otherwise.  Then  the  ijih  term  of  the  capacitance 
matrix  is  computed  by  summing  all  the  charges  on  the  jth  conductor,  i.e. 

f'tj  —  kZConductoT]  9k* 


3  Using  Preconditioned  Conjugate-Gradient 

In  order  to  solve  for  a  complete  m  x  m  capacitance  matrix,  the  n  x  n  sym¬ 
metric  matrix  of  potential  coefficients,  P,  must  be  factored  once,  usually 
into  P  =  LLt ,  and  this  requires  order  n3  operations.  Then,  as  there  are 
m  conductors,  the  factored  system  must  solved  m  times  with  m  different 
right-hand  sides,  and  this  requires  order  mn2  operations.  Since  n  is  the 
total  number  of  tiles  into  which  the  conductor  surfaces  are  cut,  m  is  neces¬ 
sarily  much  less  than  n.  Therefore,  the  n3  factorization  dominates  for  large 
problems. 

This  suggests  that  iterative  methods  might  be  more  efficient  than  direct 
factorization  for  solving  the  m  charge  distribution  problems.  In  particu¬ 
lar,  as  the  matrix  is  symmetric  and  positive  definite,  the  conjugate-gradient 
(CG)  algorithm  is  a  natural  choice[golubS3].  Unfortunately,  the  CG  algo¬ 
rithm  can  converge  slowly  when  applied  to  the  matrix  of  potential  coef¬ 
ficients,  particularly  when  the  problem  contains  widely  separated  pairs  of 
very  closely  spaced  tiles.  To  accelerate  the  convergence  of  CG,  an  attempt 
is  made  to  factor  most  of  the  part  of  the  problem  associated  with  the  closely 
spaced  tiles  directly.  To  accomplish  this,  the  smallest  cube  containing  the 
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entire  problem  is  uniformly  divided  into  a  large  number  of  cubes,  typically 
into  as  close  to  cubes  as  possible.  The  piece  of  the  potential  coefficient 
matrix  associated  with  the  tile  interactions  inside  a  cube  is  then  factored 
directly  and  used  as  a  preconditioner  to  accelerate  the  CG  algorithm.  If  the 
p  and  q  vector  in  (3)  are  reordered  so  that  tiles  contained  in  a  given  cube  are 
ordered  contiguously,  the  potential  coefficients  representing  the  interaction 
between  tiles  in  a  given  cube  will  be  blocks  on  the  diagonal  of  P.  That  is, 
p  =  P,ntracub'  +  Pintercube  where  P^tracub'  is  a  block  diagonal  matrix. 

The  CG  capacitance  extraction  algorithm  with  the  Pintracube  precondi¬ 
tioner  (PCG)  is  as  follows: 

Algorithm  1:  Preconditioned  CG  capacitance  extraction  algorithm 

Setup  Phase. 

Divide  all  the  conductors  into  a  total  of  n  tiles. 

Divide  the  tiles  into  cubes,  and  reorder  to 
make  P,„tracube  block  diagonal. 

Compute  the  Potential  Coefficient  Matrix. 
for  i  =  1  to  i  =  n 
for  j  =  1  to  j  =  n 

Compute  Pij  from  (4). 

Factor  Pintracube- 

Loop  Through  all  the  Conductors. 
for  k  =  1  to  m 

if  tile  i  is  on  conductor  k,  set  p,  =  1. 
else  pi  =  0. 

Use  PCG  to  solve  Pq  =  p. 

for  /  =  1  to  m  Ck,  =  Zktconductor,  Ik- 

Preconditioned  CG  (PCG). 

The  Setup, 
r  =  p,  q  sb  0. 

Conjugate-Gradient  Loop. 

Repeat 

Solve  Pintracube  %  —  T  • 
if  the  first  iteration  /3  =  0. 
else  /?  -  zTr/(zTr)prev. 
x  =  z  +  (3x. 
y  =  Px. 
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u  -  x^Ax- 

q  =  q  +  ax. 
t  —  t  -  ay. 
Until  Converged 


4  Acceleration  with  a  Multipole  Algorithm 

As  can  be  seen  from  examining  the  computation  in  Algorithm  1 ,  m  problems 
must  be  solved  iteratively,  and  the  major  cost  is  computing  the  matrix  P, 
and  in  each  iteration  forming  the  product  Px,  both  of  which  are  order  n2. 
This  implies  that  computing  the  capacitance  matrix  with  Algorithm  1  is 
order  mn2,  and  may  not  be  much  more  efficient  than  direct  factorization  if 
the  ratio  of  tiles  to  conductors  is  low. 

An  approach  for  reducing  the  cost  of  forming  P  and  computing  Px  in  the 
CG  algorithm  can  be  derived  by  recalling  that  if  x  is  thought  of  as  a  charge 
distribution,  Px  is  the  potential  due  to  that  charge  distribution.  To  see 
how  this  helps  simplify  the  computation  Px,  consider  two  widely  separated 
cubes,  each  with  k  tiles.  Computing  the  contributions  to  the  potentials  at 
the  center  of  each  of  the  tiles  in  the  first  cube  due  to  the  k  tile  charges  in 
the  second  cube  from  (4)  requires  k2  calculations.  If  all  the  charges  in  the 
second  cube  are  positive,  then  the  k  potential  contributions  to  the  first  cube 
can  be  computed  approximately  in  k  operations.  This  is  done  by  assuming 
the  charges  in  the  second  cube  contribute  to  potential  in  the  first  cube  like 
a  point  charge  equal  to  the  sum  of  the  charges  in  the  second  cube  located  at 
a  “center  of  mass”.  Note  that  the  accuracy  of  the  approximation  improves 
as  the  separation  between  cubes  increases. 

There  are  a  collection  of  algorithms  based  on  the  above  idea,  often  re¬ 
ferred  to  as  multipole  algorithms[rohklin86,katzenelson88,zhao87].  The  de¬ 
tails  of  the  multipole  algorithm  we  used  are  well  described  in[greengard87], 
and  only  a  very  basic  outline  will  be  given  here.  In  general,  the  potential,  il>, 
due  to  a  cube  of  point  charges  at  a  location  outside  the  radius  of  the  cube 
is  given  by  the  multipole  expansion, 

<K' ,»,*)  =  E  £  (5) 

n=Om=— n 

where  r,  9  and  <f>  are  the  spherical  coordinates  of  the  evaluation  location, 
Y™{9,<t>)  is  the  spherical  harmonic,  and  M™  is  the  multipole  coefficient, 


6 


which  can  be  computed  from  the  charge  in  the  cube  from 


(6) 

1=1 

where  p,,  a,,  and  /?,  are  the  spherical  coordinates  of  the  charge.  If 
the  evaluation  location  is  well  outside  the  cube,  then  the  potential  can  be 
accurately  computed  using  just  a  few  terms  of  the  multipole  expansion. 

Consider  a  collection  of  cubes  containing  charges  and  one  cube,  well 
separated  from  the  others,  containing  several  locations  at  which  the  potential 
must  be  evaluated.  It  is  possible  to  combine  all  the  multipole  expansions 
for  the  cubes  containing  charges  into  a  single  local  expansion  from  which 
the  potential  at  the  evaluation  points  in  the  cube  can  be  computed  quickly. 
The  local  expansion  is  given  by 


V’(r,M)  =  f;  £  L^(0,<t>)rn  (7) 

n=Om=— n 

where  r,  6  and  4>  are  the  spherical  coordinates  of  the  evaluation  location, 
an^  *  ™  are  the  local  expansion  coefficients,  which  are  computed  from  the 
combination  of  multipole  expansions  for  the  cubes  containing  charges.  Good 
accuracy  can  be  achieved  with  a  few  terms  of  the  local  expansion. 

Truncated  multipole  and  local  expansions  can  be  used  to  compute  n  po¬ 
tentials  at  n  evaluation  points  in  order  n  operations,  provided  the  charges 
and  evaluation  points  are  reasonably  separated.  To  ensure  adequate  separa¬ 
tion  and  avoid  excess  calculation,  careful  hierarchical  shifting  and  combining 
of  both  the  multipole  and  local  expansions  is  necessary,  as  is  well  described 
in  [greengard87].  With  the  computation  organized  in  this  manner,  the  mul¬ 
tipole  algorithm  can  be  used  to  compute  most  of  Px  in  Algorithm  1,  except 
the  part  due  to  interactions  between  tiles  in  a  given  cube,  and  the  tiles 
of  each  cube’s  nearest  neighbors.  This  implies  that  in  Algorithm  1,  if  the 
multipole  algorithm  is  used  to  compute  Px,  most  of  P  need  not  be  formed 
explicitly.  Note  also  that  the  part  that  must  be  computed  explicitly  includes 
P tntracube i  therefore  the  multipole  accelerated  PCG  algorithm  can  still  use 
Pintraeube  as  a  preconditioner.  Finally,  note  that  using  the  multipole  algo¬ 
rithm  to  compute  Px  implies  that  both  n2  steps  of  Algorithm  1,  forming  all 
of  P  and  computing  Px,  can  be  removed. 
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Figure  1:  Bus  Structure  Example  with  Six  Conductors 

5  Results  and  Conclusions 

The  multipole  accelerated  PCG  algorithm  was  implemented  and  tested  on 
a  simple  bus  structure  (Figure  1),  with  2,  4,  and  6  conductors.  In  Table 
1  we  report  the  total  number  of  tiles  resulting  from  the  conductor  surface 
discretization,  the  number  of  cubes  into  which  space  w as  divided,  the  time 
to  compute  capacitance  using  direct  factorization  and  PCG,  the  number  of 
iterations  to  achieve  convergence  with  PCG  and  multipole  accelerated  PCG 
(MPCG ),  and  the  relative  error  introduced  by  the  multipole  approximation. 

Much  additional  work  is  under  way  to  improve  the  efficiency  of  our 
MPCG-based  capacitance  extraction  program,  and  CPU  time  comparisons 
for  an  efficient  implementation  will  be  presented  at  the  conference.  Future 
research  includes  extending  the  approach  to  piecewise-constant  dielectrics 
and  problems  with  ground  planes. 

The  authors  would  like  to  thank  David  Ling  and  Albert  Ruehli  of  the 
I.B.M.  T.  J.  Watson  Research  Center  for  the  many  discussions  that  led  to  the 
approach  presented  here,  as  well  as  their  help  along  the  way.  In  addition  we 
would  like  to  acknowledge  the  helpful  discussions  with  Jacob  Katzenelson, 
and  finally  we  thank  the  many  members  of  the  MIT  Custom  Integrated 
Circuits  group  for  their  help  and  encouragement. 
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2  Cond. 

4  Cond. 

6  Cond. 

tiles 

216 

720 

1512 

cubes 

64 

64 

64 

direct  time 

67 

2520 

20160 

PCG  time 

65 

653 

2613 

PCG  iters 

T 

1 

10 

12 

MPCG  iters 

8 

11 

14 

MPCG  rel.  err. 

0.002 

0.001 

0.002 

Table  1:  Comparison  of  Extraction  Methods 

Agency  contract  N00014-87-K-825,  and  grants  from  IBM  and  Analog  De¬ 
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